Notebook para criar a relação periodo luminosidade para as RRAB do ogle III
O arquivo de possui os seguintes dados:
StarID Imédio Vmédio PeriodoOGLE RA DE PeriodoCE
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
#load file
data = np.loadtxt('/home/glauffer/Dropbox/FURG/final_project/data/ID_mag_Period_RA_DE.dat', dtype = float, usecols = (1,3,6))
#data = np.genfromtxt('/home/glauffer/Dropbox/FURG/final_project/data/ID_mag_Period_RA_DE.dat', dtype=str)
#mean_I = data.T[1]
#mean_I
plt.clf()
#plt.plot(np.log(data.T[2]), data.T[0], 'bo')
plt.plot(data.T[2], data.T[0], 'bo')
plt.gca().invert_yaxis()
plt.xlim(0.45,0.7)
plt.ylim(21, 18)
Out[1]:
o modulo de distancia, "Distance Modulus", para a LMC, de acordo com Schaefer (2008) possui valor $\mu \approx 18.5 \pm 0.1$ mag
A relação periodo luminosidade (citar a henrieta) é uma relação do tipo $$ m_i = a \log P_i + b $$
Para melhor fitar uma curva deste tipo nos dados, vou utilizar uma regressão de mínimos quadrados para achar $a$ e $b$ que melhor fitam a curva aos dados.
In [2]:
A = np.column_stack((np.log(data.T[2]), np.ones_like(data.T[2])))
mags = data.T[0]
x, res, rank, s = np.linalg.lstsq(A,mags)
print('m = ', x[0], '*log P + ', x[1])
In [15]:
plt.clf()
plt.plot(A.T[0], mags, 'bo', label = 'data', ms=0.2)
plt.plot(A.T[0], x[0]*A.T[0]+x[1], '-r', label = 'fit')
plt.plot(A.T[0], (-1.62)*A.T[0]+18.391, '-k', label = 'ogle3')
plt.legend()
plt.xlabel(r"$\log P$")
plt.ylabel(r"$<I> mag$")
plt.gca().invert_yaxis()
plt.xlim(np.log(0.45),np.log(0.7))
plt.ylim(21, 18)
#plt.xlim(-3, 1)
#plt.savefig('pl_relation_RRAB.png', dpi = 100)
Out[15]:
In [4]:
#obtendo m_fit
m_fit = x[0]*A.T[0] + x[1]
residual = mags - m_fit
residual
plt.plot(A.T[0], residual, 'ro', label = 'residual')
plt.legend()
plt.gca().invert_yaxis
Out[4]:
ver sobre peso nos residuos
In [5]:
#utilizando os residuais para determinar distancia (nao sei se esta correto)
d = 10**(0.2*(residual+5))
d.mean()
#calculando o modulo de distancia pela diferença entre o artigo e o fit (tambem nao sei se esta certo)
artigo = (-1.453)*A.T[0]+14.828
res1 = artigo - m_fit
mu = 14.828 - x[1]
res1
Out[5]:
In [5]: